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Abstract 

In  previous  work  we  have  suggested  obtaining  rational  interpolants  of  a function  / by 
attaching  optimally  placed  poles  to  its  interpolating  polynomials.  For  a large  number  of 
interpolation  points  these  polynomials  are  well-known  to  be  good  approximants  only  if 
the  nodes  tend  to  cluster  near  the  endpoints  of  the  interval,  as  with  Cebysev  or  Legendre 
points.  In  practice,  however,  one  would  prefer  to  have  them  closer  to  equidistant.  This 
will  in  particular  be  the  case  when  the  difficult  portion  of  / lies  well  within  the  interior 
of  the  interval,  or  when  approximating  derivatives  of  /,  as  in  the  solution  of  differential 
equations.  To  address  this  difficulty,  we  use  here  a conformal  change  of  variable  to  shift 
the  points  from  the  Cebysev  position  toward  a more  equidistant  distribution  in  a way 
that  should  maintain  the  exponential  convergence  when  / is  analytic.  Numerical  examples 
demonstrate  the  resulting  improvement  in  the  quality  of  the  approximation. 


1 Introduction 

We  are  concerned  here  with  rational  approximation  of  a continuous  function  / on  an 
interval  [a,  b ],  which  we  may  take  as  [—1, 1]  =:  I,  after  a linear  change  of  variable  when 
necessary.  We  further  assume  that  the  approximant.  r should  interpolate  / between  a 
finite  number,  say  N + 1,  of  distinct  points  (nodes)  x0,  aq , . . . , aqv  in  I.  In  a similar  way 
as  in  [5],  r will  be  constructed  by  attaching  a certain  number  of  poles  to  an  interpolating 
polynomial. 

In  some  applications,  such  as  the  numerical  solution  of  two-point,  boundary  value 
problems  (see,  e.g.,  [6]),  one  may  choose  the  points  more  or  less  at  will;  in  that  case, 
one  will  place  them  so  as  to  reach  the  best  compromise  between  two  often  conflicting 
goals:  points  good  for  interpolation,  on  one  side,  and  points  favourable  for  the  condition 
of  the  problem  to  be  solved,  on  the  other  side.  In  [5],  we  have  considered  equidistant 
and  Cebysev  points,  the  first  for  their  regularity,  the  second  for  the  condition  of  the 
interpolation  and  for  the  fast  convergence  of  the  interpolant  for  very  smooth  functions. 
For  the  solution  of  two-point  boundary  problems  in  [6]  we  have  merely  used  Cebysev 
points. 
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There  is  in  general  no  reason  besides  the  problem  condition  for  accumulating  the 
nodes  toward  the  boundary,  as  with  Cebysev  or  Legendre  points.  Moreover,  one  of  the 
reasons  for  using  rational  instead  of  polynomial  interpolation  is  its  better  suitability  for 
approximating  functions  with  large  slopes.  Here  too,  shifting  the  points  away  from  the 
center  may  not  be  appropriate. 

Another  odd  consequence  of  accumulating  interpolation  points  toward  the  extremities 
is  the  consequent  ill-conditioning  of  the  derivatives  of  the  interpolating  polynomials  [7, 1] . 

This  worsens  the  stability  properties  of  time-stepping  in  the  solution  of  time  evolution 
problems  with  the  method  of  lines  [13]  as  well  as  the  convergence  of  iterative  methods 
for  solving  discretized  stationary  problems  [3], 

To  address  these  difficulties,  we  will  take  advantage  here  of  the  fact  that  the  fast  con- 
vergence of  the  interpolant  can  be  maintained  while  shifting  the  points  with  a conformal 
map  g (independent  of  N ) toward  an  equidistant  position.  This,  however,  requires  an 
important  change  to  the  method  in  [5],  because  this  point  shift  ruins  the  exponential 
convergence  of  the  Cebysev  interpolating  polynomial.  We  therefore  use  here  as  the  start- 
ing interpolant  the  polynomial  interpolating  /(#-1)  in  the  domain  of  the  inverse  #-1  of 
the  conformal  map  employed  for  the  point  shift,  and  attach  poles  to  this  polynomial. 

Section  2 reviews  the  formulae  and  advantages  of  shifting  Cebysev  points  conformally 
toward  the  center  of  the  interval  when  interpolating  functions,  and  Section  3 briefly  re- 
calls the  method  of  optimally  attaching  poles  to  the  interpolating  polynomial  introduced 
in  our  earlier  work.  In  Section  4 we  describe  how  to  take  advantage  of  the  better  condi- 
tioning of  derivatives  induced  by  the  conformal  point  shift;  the  corresponding  practical 
improvements  are  finally  documented  with  numerical  examples. 


2 Rational  interpolation  with  a variable  change  for  point  shifts 

Let  Vm  and  respectively,  denote  the  linear  space  of  all  polynomials  of  degree  at 

most  m and  the  set  of  all  rational  functions  with  numerator  in  Vm  and  denominator  in 
Vn\  furthermore,  denote  by  fk  the  interpolated  values  f(xk ),  k = 0(1) N,  of  /.  Then,  the 
unique  polynomial  p £ Vn  that  interpolates  / at  the  x^s, 

N 

p(x)  = ^2  fkLk{x ),  Lk(x)  :=  PJ(x  - Xi)  I JJ(xfc  - xf), 

k— 0 i^k  / iy£k 


can  be  written  in  its  barycentric  form  [9] 


N 

p(x)  = Y2 


k= 0 


Wk 

x-xk 


Wk 

X Xk 


(2.1) 


where  the  so-called  weight  wk  corresponding  to  the  point  xk  is  given  by 


N 


Wk  :=  1 / JJ  (xk-Xi). 

/ i=0,  i^k 


Despite  its  appearance,  (2.1)  determines  a polynomial  of  degree  at  most  N:  the  wk 
are  precisely  the  numbers  which  guarantee  this  [4].  By  choosing  other  wks,  a rational 
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interpolant  is  constructed. 

The  barycentric  formula  has  several  advantages  over  other  representations  of  the 
interpolating  polynomial  ([4]  p.  357).  One  of  them  is  the  fact  that  the  weights  appear  in 
both  the  numerator  and  the  denominator,  so  that  they  can  be  divided  by  any  common 
factor.  For  example,  simplified  weights  for  Cebysev  points  of  the  first  kind  cos(j>k, 

where  (j>k  27^+1) ^ and  ^ = 0 are  given  by  = (-l)fc sin^.  ([9]  p.  249), 

while  for  the  Cebysev  points  of  the  second  kind  xk  ’ :=  cos  k - which  will  be  used  here 
- one  simply  has  Salzer’s  formula  ([9]  p.  252) 


which  will  be  used  here 


= (-!)*%, 


-(f 


k = 0 or  k = N, 
otherwise. 


These  points  are,  together  with  Legendre’s,  the  most  used  nodes  for  global  polynomial 
interpolation  and  large  N.  They  achieve  exponential  convergence  of  p toward  / if  the 
latter  is  analytic  in  an  ellipse  Ep  with  foci  at  ±1  and  sum  of  its  axes  equal  to  2p,  p > 1. 
However,  this  fast  convergence  comes  at  the  cost  of  a concentration  of  the  nodes  in 
the  vicinity  of  the  extremities  of  I.  As  mentioned  above,  this  accumulation  may  have 
drawbacks,  such  as  poor  spreading  of  the  information  about  / over  the  interval  and 
ill-conditioning  of  the  derivatives  near  the  endpoints. 

With  a suitable  choice  of  the  interpolant,  one  may  conformally  shift  the  nodes  to- 
ward an  equidistant  position  (though  not  all  the  way)  without  losing  the  exponential 
convergence.  For  that  purpose,  one  considers,  beside  the  x-space  in  which  / is  to  be 
approximated,  another  space,  denoted  by  y,  say,  and  the  N + 1 Cebysev  points  of  the 
second  kind 

(2) 

Vk  = 4 

in  the  interval  J :=  [-1,1]  in  this  y- space.  Let  g be  a conformal  map  from  a domain  V\ 
containing  J (in  the  y- space)  to  a domain  V2  containing  I (in  the  x-space);  moreover, 
suppose  that  / is  a function  X>2|— *•  C such  that  the  composition  fog:  C is  analytic 

in  an  ellipse  Ep,  as  defined  above.  With  this  map  we  may  define  new  interpolation  points 
on  I,  Xk  = g{yk),  as  well  as  the  conformal  transplantation  F(y)  :=  f(x)  [10]  of  / into 
the  y- space. 

Then,  with  the  polynomial  interpolating  F(y)  at  the  yk 
N N 

An{v)  :=^2F{yk)Lk{y)  = '}Tf(xk)Lk{g~1{x))  -■  aN(x),  (2.2) 


one  has 


lajvOc)  - f(x)\  = 0(p  N),  x € [—1, 1]. 


Rational  interpolation  with  all  poles  prescribed  is  very  simple  in  the  barycentric  setting 
[5]:  the  P poles  z;  are  attached  to  (2.1)  by  replacing  Wk  with 


bk  = wkdk , dk  :=  JJ(xfc  - zt ). 
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If  N > P this  results  in  a rational  interpolant  in  7Zn,p  with  poles  at  Zi,  i = 

(when  such  an  interpolant  exists,  see  [5]). 

Remark  2.1  Exponential  convergence  of  interpolation  at  the  shifted  points  is  also  at- 
tained with  the  rational  function  given  by  (2.1)  with  Wk  = [2].  However,  this  is 

in  general  a rational  function  in  Rn,v,  v > N — P:  there  is  not  enough  defect  in  the 
denominator  degree  for  the  weights  wKk  'dk  to  warrant  the  presence  of  the  P poles  Zi. 

We  then  use  as  the  starting  interpolant  to  which  we  attach  the  poles  m in  the 
y- space.  This  yields 


N 


Wk 


II<*  - Vi) 


N 


i(xk)-9  Hzi)) 


R(y)  := 


V-'  i= 1 , \T'  i= 1 

^ V~Vk  k (rfy  g-1(x)-g-1(xk) 


-fk 


k= 0 


N 


Wk  - Vi) 


WkU(9  1(xk)-g  1(zi)) 


r(x). 


N 


E 

k= 0 


j=l 


y-yk 


2 9 1{x)-9  1{xk) 


If  a rational  interpolant  with  these  poles  exists,  it  is  given  in  the  y-space  by  R,  and  r is 
a rational  function  in  the  argument  g~1(x).  Its  poles  are  at  z,  — g{vi). 


3 Construction  of  the  optimal  interpolant 

Our  method  consists  in  optimizing  the  position  of  the  vfs  so  as  to  minimize 

IIR-ClooHk-ZIU 

as  described  in  §3  of  [5].  Optimal  vf  s always  exist,  but  these  are  not  unique  in  general. 
Whether  the  optimal  R is  unique  is  an  open  question;  however,  for  every  optimized  pole 
Vi  an  indicator  may  be  calculated  which,  if  nonzero,  guarantees  that  ty  is  indeed  a pole 
of  R. 

In  the  practical  computations  documented  in  §5  the  optimization  of  the  iVs  was 
performed  using  the  same  two  algorithms  as  in  [5]:  for  small  N we  used  a discrete  differ- 
ential correction  algorithm  according  to  [11],  while  for  larger  N the  simulated  annealing 
method  of  [8]  was  applied.  Both  methods  will  in  principle  locate  a desired  global  maxi- 
mum. The  first  method  achieves  it  in  a systematic  and  guaranteed  way  evaluating  the 
error  not  continuously  but  on  a fine  grid;  the  simulated  annealing  method  cannot  be 
guaranteed  to  find  the  global  extremum  but,  when  used  for  an  extensive  search,  will 
produce  a reasonable  approximation  of  it. 

As  mentioned  in  [5] , our  way  of  attaching  poles  to  the  interpolating  polynomial  has  a 
very  nice  property:  the  approximation  error  can  only  decrease  or  at  worst  stay  Constant 
with  a growing  number  of  poles,  this  in  sharp  contrast  with  classical  rational  interpola- 
tion; when  a new  unknown,  say  Vj,  is  added  to  the  set  of  variables,  {iq, . . . , v j _ i } , the 
optimal  values  of  the  latter  are  a feasible  vector  for  the  higher  dimensional  optimization. 

Let  us  conclude  this  section  with  a comment  on  the  use  of  the  nomenclature  “at- 
taching the  poles”.  In  classical  rational  interpolation,  the  poles  of  the  interpolant  are 


Jean-Paul  Beirut  and  Hans  D.  Mittelmann 


determined  by  the  data.  There  too,  however,  one  sometimes  wishes  to  prescribe  the  loc- 
ation of  the  poles  (with  corresponding  decrease  of  the  number  of  degrees  of  freedom): 
many  authors  then  speak  of  “assigning”,  or  “prescribing”  the  poles.  In  that  sense  one 
cannot  “assign”  poles  to  a polynomial,  which  obviously  cannot  have  poles.  We  thus 
start  with  the  interpolating  polynomial  and  its  poles  at  infinity  and  make  it  a rational 
interpolant  by  bringing  the  poles  into  an  optimal  position  in  C.  We  call  this  procedure 
“attaching  poles”,  to  distinguish  it  from  the  process  of  forcing  a rational  function  to 
have  a pole  at  a particular  place. 

4 Derivatives  of  the  optimal  interpolant  with  shifted  points 

As  mentioned  in  §1,  one  of  the  reasons  for  shifting  the  points  from  their  Cebysev  position 
toward  the  interior  of  the  interval  is  the  improvement  of  the  condition  of  the  derivatives 
resulting  from  such  a shift.  Besides  r,  we  will  evaluate  also  r'  and  r"  as  approximants  of 
/',  resp.  /",  and  estimate  ||r  - /'||oo  and  ||r  - /"||oo- 

Schneider  and  Werner  [14]  have  noticed  that  every  rational  interpolant  R £ R-n,n, 
written  in  its  barycentric  form 


y — 2/a- 


can  easily  be  differentiated.  The  formulae  for  the  first  two  derivatives  read 


R'(y)  = 


N IN 

E— «[»•»*] /E— , 
; t0y-y*  /; toy~yk 

-(  E uk  R[yi,  2/a-])  j ui , 


Vi,  i = 0{1)N, 


y = yi 


2E^[^WE 


R"(y)  = 


,y-yk 


y-yk 


y^Vi,  i = 0(1) A, 


-2 ( £ v-kR[vuVi, yd) / «»!  y = yu 


ith  R[z, z, j/fc]  = - ■ The  chain  rule  then  yields,  for  r(x)  = R[g  1(x)), 

/(x)  = R'(y)  ■ [»-■«]'  = = |57^5pR"(»>  - (»)•  (4.1) 


Specifically,  in  our  calculations  we  have  used  the  map  suggested  by  Kosloff  and  Tal-Ezer 

[12], 

, . arcsin(ay) 

g(y)  = r-^,  0 < a < 1. 

arcsin  a 
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a 

P = 0 

P = 2 

warn 

P = 6 

P = 8 

0.0 

6.37e  - 5 

1.42e  - 6 

5.83e  - 8 

9.38e  - 9 

1.30e  - 9 

0.5 

3. lie  — 5 

6.69e  - 7 

2.48e  - 8 

4.21e  - 9 

4.23e  - 10 

wmnm 

8.06e  — 6 

1.60e  - 7 

5.50e  - 9 

9.47e  - 10 

1.27e  - 10 

0.9 

1.12e  - 6 

1.97e  - 8 

■Xli'LTi'lffpqgfyp 

2.05e-  11 

0.95 

2.78e  - 7 

4.47e  - 9 

1.29e  - 10 

1.36e  - 11 

3.82e  - 12 

0.96 

1.85e  - 7 

2.93e  - 9 

8.27e-  11 

4.20e  - 12 

3.88e  - 12 

Tab.  1.  Errors  when  approximating  / with  increasing  P and  a in  Example  1. 

In  the  limiting  cases,  a — > 0 keeps  the  points  at  their  Cebysev  position,  whereas  a — > 1 
renders  them  equidistant.  The  derivatives  of  g are  given  by 


g'{y)  = 


so  that  in  (4.1) 


- {ay)2  ’ 

9"{V) 

[g'{y)f 


- (arcsin2  a)y. 


Numerical  evidence 


We  now  report  on  practical  computations,  performed  on  two  examples,  which  demon- 
strate the  efficiency  of  point  shifts  for  improving  the  rational  interpolants  with  optimized 
denominators.  These  examples  share  the  property  that  the  difficult  part  of  / lies  in  the 
center  of  I , so  that  the  shift  of  the  points  toward  a more  equidistant  position  naturally 
improves  the  quality  of  the  information  provided  to  the  interpolation  method. 


P = 0 

P = 2 

5.27e  - 3 

1.26e  - 3 

2.67e  - 3 

5.87e  - 4 

7.47e  - 4 

1.49e  - 5 

1.14e  — 4 

2. Ole  - 6 

2.97e  - 5 

4.99e  - 7 

2. Ole  — 5 

3.24e  - 7 

= 4 


4.85e  - 6 


P = 6 


8.69e  - 7 


P = 

8 

1.40e 

-7 

4.63e 

-8 

1.30e 

-8 

2.16e 

-9 

4.52e  - 

- 10 

4.70e  - 

- 10 

Tab.  2.  Errors  when  approximating  /'  with  increasing  P and  a in  Example  1. 

The  sup-norm  ||  ||oo  has  thereby  been  estimated  by  considering  the  1000  equally 
spaced  points  xe  = — | + f^xi  ^ = 1(1)1000,  on  the  interval  [—5/4, 5/4]  and  computing 
the  maximal  absolute  value  of  the  error  at  those  xe  lying  in  [—1, 1]. 

Example  5.1  We  have  first  revisited  Example  3 of  [5],  which  displays  in  the  center  of 
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I a slope  increasing  with  a positive  parameter,  here  denoted  by  e, 


/( x)  = COS  7TX  + 


erf(5.x) 
erf(<5)  ’ 


6 — vC5e, 


where  erf  denotes  the  error  function  (see  [5]  for  a graph). 


In  Table  1 we  give  the  results  obtained  with  e = 500  and  N = 81.  increasing  numbers 
P of  poles  and  increasing  a.  Tables  2 and  3 display  the  same  information  for  the  appro- 
ximation of  f and  f"  with  r'  and  r"  as  given  by  the  formulae  (4.1).  The  combination 
of  extra  poles  and  a point  shift  brings  about  7 digits  of  accuracy,  where  the  point  shift 
alone  makes  only  for  2-3.  The  improvement  in  the  derivatives  is  especially  remarkable: 
the  error  in  the  second  derivative  decreases  from  the  useless  value  of  9.26  to  about  10~7! 


a 

P = 0 

<N 

II 

ft, 

■Bl 

wmmm 

P = 6 

0.0 

9.26 

4.05e  - 2 

4.82e  - 4 

7.85e  - 5 

1.46e  - 5 

0.5 

4.26 

2.07e  - 2 

2.18e  - 4 

3.75e  - 5 

4.91e  - 6 

0.75 

9.50e  - 

1 

5.48e  - 3 

6.25e  - 5 

9.53e  - 6 

1.26e  - 6 

0.9 

9.30e  - 

2 

6.49e  — 4 

8.86e  - 6 

4.93e  - 7 

2.34e  - 7 

0.95 

1.59e  - 

2 

1.23e  — 4 

1.88e  - 6 

1.75e  - 7 

5.31e  - 8 

0.96 

9.18e  - 

3 

7.36e  - 5 

1.29e  - 6 

6.00c  - 8 

5.57e  - 8 

Tab.  3.  Errors  when  approximating  /"  with  increasing  P and  a in  Example  1. 


Example  5.2  Example  3 in  [5]  has  demonstrated  that  the  attachment  of  poles  may  be 
very  effective  in  improving  the  approximation  of  oscillatory  functions.  Here  we  change 
the  function  to 

h(x)  = e~ax  sin6x,  a > 0,  b > 0, 
so  that  the  most  oscillatory  part  lies  in  the  center  of  the  interval. 

Results  with  a — 5,  b = 25,  N = 31,  P = 0 and  P — 2 are  given  in  Table  4.  In  contrast 
with  the  preceding  example,  here  the  point  shift  brings  much  more  improvement  than  the 
attachment  of  poles,  about  6-7  digits,  an  especially  heartening  fact  for  the  derivatives, 
to  which  the  interpolants  without  shift  are  useless  approximants. 

Acknowledgement:  The  authors  wish  to  thank  Peter  Graves-Morris  for  his  comments 
which  have  enhanced  the  present  text. 
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a 

h 

h" 

P = 0 

P=:  2 

P = 0 

P = 2 

P = 0 
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1.66e  - 2 

8.68e  - 4 
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Tab.  4.  Change  in  the  errors  induced  by  the  introduction  of  two  poles  in  Example  2. 
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